knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',
                      echo = TRUE, message = FALSE, warning = FALSE)

library(raster)
library(data.table)
library(sf)


source('https://raw.githubusercontent.com/oharac/src/master/R/common.R')  ###
  ### includes library(tidyverse); library(stringr); dir_M points to ohi directory

dir_git <- '~/github/spp_risk_dists'

### goal specific folders and info
dir_data  <- file.path(dir_git, 'data')
dir_o_anx <- file.path(dir_O, 'git-annex/spp_risk_dists')

source(file.path(dir_git, 'setup/common_fxns.R'))

### Gall-Peters doesn't have an EPSG?
gp_proj4 <- '+proj=cea +lon_0=0 +lat_ts=45 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs'

1 Summary

Create a set of maps of the distribution of biodiversity intactness - all species assessed and mapped by IUCN. These maps are generated at 10 km² resolution in a Gall-Peters projection. These maps will be generated using all available species:

Mean risk
Variance of risk
Number of species for mean/var calculations
Number of species categorized as “threatened” (i.e. VU, EN, CR)
Mean trend
Number of species for trend calculations

A selection of these maps will be generated for taxonomic groups and range sizes in a separate Rmd.

Future iterations may include:

Range-rarity-weighted mean and variance of risk
Range rarity-weighted species richness

2 Data Sources

IUCN Red List spatial data download IUCN Red List API Gina Ralph (IUCN)

3 Methods

3.1 Spatial distribution of risk: comprehensively assessed spp

From the 1a and 1b biodiversity risk map scripts, gather the rasters for the various maps:

Mean risk (unweighted and range-rarity-weighted)
Pct threatened (unweighted and range-rarity-weighted)

risk_un_rast    <- raster(file.path(dir_git, 'output', 'mean_risk_raster_comp.tif'))
risk_rr_rast    <- raster(file.path(dir_git, 'output', 'mean_rr_risk_raster_comp.tif'))
threat_un_rast  <- raster(file.path(dir_git, 'output', 'pct_threat_raster_comp.tif'))
threat_rr_rast  <- raster(file.path(dir_git, 'output', 'sr_rr_pct_threat_raster_comp.tif'))
n_spp_rast <- raster(file.path(dir_git, 'output', 'n_spp_risk_raster_comp.tif'))
n_rr_rast  <- raster(file.path(dir_git, 'output', 'sr_rr_risk_raster_comp.tif'))

3.1.1 And now, the maps

four_panel <- function(map1_rast, map2_rast,
                       limits = c(0, 1),
                       colors, values,
                       labels, breaks,
                       plot_labs = c('A', 'B')) {

  land_poly <- sf::read_sf(file.path(dir_git, 'spatial/ne_10m_land/ne_10m_land.shp')) %>%
    st_transform(gp_proj4)
  
  map_df <- data.frame(val_1 = values(map1_rast),
                       val_2 = values(map2_rast)) %>%
    cbind(coordinates(map1_rast)) %>% 
    filter(!is.na(val_1))
  
  map1 <- ggplot(map_df) +
    geom_raster(aes(x, y, fill = val_1), show.legend = FALSE) +
    geom_sf(data = land_poly, aes(geometry = geometry), 
            fill = 'grey96', color = 'grey40', size = .10) +
    ggtheme_map() +
    theme(plot.margin = unit(c(.2, 0, .1, .5), units = 'cm')) +
    coord_sf(datum = NA) + ### ditch graticules
    scale_fill_gradientn(colors = colors, values = values, limits = limits,
                         labels = labels, breaks = breaks,
                         guide  = guide_colourbar(label.position = 'left',
                                                  label.hjust = 1)) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0))
  
  map2 <- ggplot(map_df) +
    geom_raster(aes(x, y, fill = val_2), show.legend = FALSE) +
    geom_sf(data = land_poly, aes(geometry = geometry), 
            fill = 'grey96', color = 'grey40', size = .10) +
    ggtheme_map() +
    theme(plot.margin = unit(c(.1, 0, .2, .5), units = 'cm')) +
    coord_sf(datum = NA) + ### ditch graticules
    scale_fill_gradientn(colors = colors, values = values, limits = limits,
                         labels = labels, breaks = breaks,
                         guide  = guide_colourbar(label.position = 'left',
                                                  label.hjust = 1)) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0))

  ###################
  
  ### set up a dataframe of values to craft a color bar using geom_segment
  colorbar_df <- data.frame(x = seq(0, 1, .001), y = -1)

  dens1 <- ggplot(map_df) +
    ggtheme_plot() +
    geom_vline(xintercept = mean(map_df$val_1, na.rm = TRUE)) +
    geom_segment(data = colorbar_df, 
                 aes(x = x, xend = x, color = x), 
                 y = 0, yend = -10, size = 1,
                 show.legend = FALSE) +
    geom_density(aes(x = val_1, ..scaled..), alpha = .3, size = .25, fill = 'grey30') +
    scale_color_gradientn(colors = colors, values = values, limits = limits,
                          labels = labels, breaks = breaks) +
    theme(axis.text.x  = element_blank(), 
          axis.title = element_blank(),
          panel.grid.major.x = element_blank(),
          plot.margin = unit(c(1, 0, 1, 0), units = 'cm')) +
    scale_x_continuous(labels = labels, breaks = breaks, limits = limits,
                       expand = c(0, 0)) +
    scale_y_continuous(limits = c(-.4, 1)) +
    coord_flip()
    

  dens2 <- ggplot(map_df) +
    ggtheme_plot() +
    geom_vline(xintercept = mean(map_df$val_2, na.rm = TRUE)) +
    geom_segment(data = colorbar_df, 
                 aes(x = x, xend = x, color = x), 
                 y = 0, yend = -10, size = 1,
                 show.legend = FALSE) +
    geom_density(aes(x = val_2, ..scaled..), alpha = .3, size = .25, fill = 'grey30') +
    scale_color_gradientn(colors = colors, values = values, limits = limits,
                          labels = labels, breaks = breaks) +
    theme(axis.text.x  = element_blank(), 
          axis.title = element_blank(),
          panel.grid.major.x = element_blank(),
          plot.margin = unit(c(1, 0, 1, 0), units = 'cm')) +
    scale_x_continuous(labels = labels, breaks = breaks, limits = limits,
                       expand = c(0, 0)) +
    scale_y_continuous(limits = c(-.4, 1)) +
    coord_flip()
 
  panel_top <- cowplot::plot_grid(map1, dens1, 
                                  axis = 'b',
                                  rel_widths = c(5, 1))
  panel_btm <- cowplot::plot_grid(map2, dens2, 
                                  axis = 'b',
                                  rel_widths = c(5, 1))
  
  four_panel <- cowplot::plot_grid(panel_top, panel_btm, 
                                   labels = plot_labs,
                                   nrow = 2, ncol = 1, align = 'v')

}

### aggregate rasters for faster testing:
# risk_un_rast <- risk_un_rast %>% aggregate(10)
# risk_rr_rast <- risk_rr_rast %>% aggregate(10)

mean_four_panel <- four_panel(risk_un_rast, risk_rr_rast,
                              colors = risk_cols,
                              values = risk_vals,
                              labels = risk_lbls,
                              breaks = risk_brks,
                              plot_labs = c('A', 'B'))

### aggregate rasters for faster testing:
# threat_un_rast <- threat_un_rast %>% aggregate(10)
# threat_rr_rast <- threat_rr_rast %>% aggregate(10)

threat_four_panel <- four_panel(threat_un_rast, threat_rr_rast,
                                colors = thr_cols,
                                values = thr_vals, 
                                breaks = thr_brks,
                                labels = thr_lbls,
                                plot_labs = c('C', 'D'))

four_maps <- cowplot::plot_grid(mean_four_panel, threat_four_panel, 
                                nrow = 1, ncol = 2, align = 'v')

fig1_path <- file.path(dir_git, 'ms_figures/fig1_comp_assessed.png')
ggsave(plot = four_maps,
       filename = fig1_path,
       width = 10, height = 5, units = 'in', dpi = 300)

3.1.2 Fig 1: Biodiversity risk and threatened species: comprehensively assessed spp

3.2 Spatial distribution of risk: ALL spp